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Administrative Notes: 


The grant funds were released only on November 1,1995 with back dating of the grant to 
August 1,1995. A new graduate student (Mr. Luan Nguyen) could be appointed on the grant 
only on January 2, 1996. The work reported here was done by the PI and a former graduate 
student. It will be continued by the PI and the new graduate student for the remaining duration of 
the project 

Technical Progress - Summary : 

Work Completed: 

In our earlier work, we developed mathematical models that can describe the behavior of 
free falling wavy films at high Reynolds numbers. One major objective of our proposed work 
was to extend these models to include the co or counter-current flow of gas. In the first phase of 
our work, we have completed the extension and analysis of the model for counter-current gas 
flow. The details of the model and the analysis are given in the Appendix to this report. We are 
presently in the process of writing up the results of this study as a publication in a refereed 
journal. 

Work to be Completed: 

The PI and the new graduate student are now working on the extension of the boundary 
layer models to microgravity conditions. We hope to complete this extension as well as obtain 
some theoretical results that may be used to compare with experimental data. We also plan to 
complete the study of heat (and mass) transfer enhancement studies in wavy films at high 
Reynolds numbers. The results of these two extensions will be reported in the final report. 



Semi-annual Progress Report: 


Grant # NAG 9-854: 

Heat and Momentum Transfer Studies in High Reynolds Number Wavy Films at Normal and 

Reduced Gravity Conditions 


1. Introduction: 

The most dramatic effect of the simultaneous flow of gas and liquid in pipes is the greatly 
increased transport rates of heat, mass and momentum. In practical situations this enhancement 
can be a benefit or it can result in serious operational problems. For example, gas-liquid flow 
always results in substantially higher pressure drop and this is usually undesirable. However, 
much higher heat transfer coefficients can be expected and this can obviously be of benefit for 
purposes of design. Unfortunately, designers know so little of the behavior of such two phase 
systems and as a result these advantages are not utilized. In the first phase of our work, we 
examined the effect of the gas flow on the liquid film when the gas flows in the countercurrent 
direction in a vertical pipe at normal gravity conditions. 

Only a few theoretical studies dealt with wavy motion with countercurrent or cocurrent 
gas flow. Kapitza [1] studied the effect of countercurrent or cocurrent gas flow in the momentum 
equation for the interfacial shear stress, which arises from the gas-liquid interaction. The most 
significant contribution to wave modeling for horizontal gas-liquid flow in thin film has been 
presented by Miya et al.[2], A solitary wave was considered and a shock condition at the wave 
front was utilized. This model satisfactorily predicts the interfacial shape, wall shear stress, 
pressure drop, and other wave variables. The approach is still rather approximate since 
experimental data is needed as an input to provide a solution. Brauner et al.[3] used an integral 
approach to model waves observed in cocurrent downward annular gas-liquid flow. 

The experimental studies of countercurrent flow have been made by many researchers. 
Important parameters for the description of countercurrent flow are pressure drop, wall shear 
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stress and interfacial wave structure. Pressure drop measurements have been made by Feind[4] , 
Dukler and Smith [5], and Zabaras[6]. For adiabatic flow, the pressure drop slowly increases 
before the flooding point as the gas rate is increased under constant liquid rate. Most studies of 
interfacial wave structure have been restricted in measurements of the time average film 
thickness (Feind[4], Dukler and Smith[5]). The results of these investigators show that the mean 
film thickness is not significantly affected by the upward flow of the gas, although no 
quantitative expressions were presented. Experimental results indicate that the mean film 
thickness with countercurrent gas flow increases by approximately 10% to 20% over the film 
thicknesses with zero gas flow. Wave frequencies and wave velocities in vertical countercurrent 
flow are given by Hewitt and Wallis [7] and Zabaras [6]. These investigators report that there is 
little influence of the countercurrent air flow on the wave velocities in falling films. 

2. Model Formulation 

We consider a viscous fluid film flowing down a vertical plate under gravity with a 
countercurrent stream of gas phase adjoining the free surface. The particular flow situation 
considered is one for which there is no exchange of droplets between the film and gas. We 
assume that the film thickness is considerably smaller than any length scale in the mean flow 
direction and all the variables are independent of the transverse coordinate. The two-dimensional 
Navier-Stokes equations for an incompressible fluid are 

u*Vu = - — Vp + g + vV 2 u , (2.1) 

at p 

Vu =0, (2.2) 

where u = (u,v) is the velocity field, and g = (g, 0). The coordinate system is chosen such that the 
y direction is normal to the wall and x direction is aligned with gravity (Figure 1). The pressure 
in the gas phase is p g , the kinematic and interface stress conditions can be written as 

v(h) = u(h) ^ @ y = h(x,t) , 


(2.3) 
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(|r + s)( 1 - h < 2 ) + 2 (|- s) h * = b ® y = h(x ' ,) • 


dv du' 


/du 3v^h x (Ai n x o dv l 

p-p 8 + ^ + 5TJi^2- 2 ^ - 2m ¥T^7' 

h xx 


du h* 


dv 1 


(l+hx 2 ) 3/2 


= 0 @ y = h , 


where 



and Xj is the interfacial shear stress. 

There is also the no-slip condition at the wall, 

u = 0 @ y = 0 . 


(2.4) 


(2.5) 


( 2 . 6 ) 


For the gas phase, a shallow fluid assumption is made to simply the relation between p g and film 
thickness h. The justification for this assumption for the gas flow is not so good as for the liquid 
film flow since the thickness of gas layer is much greater than the thickness of the liquid layer. 
Support for this assumption can be found elsewhere [2]. 



By using the shallow gas assumption, a momentum balance and mass balance can be developed 
for the gas phase. For the case of a countercurrent stream of gas phase adjoining the free surface, 
we have 




dx 


2tj dp t 

D-2h ' dx 


(2.7) 


and 


(D-2h)(U a +V w ) = (D-2h N )(U ao 4u Ni ) , (2.8) 

where D is the diameter of the column, p g , the gas density, U a , the local average gas velocity and 
U ao , the average gas velocity over the base film, V w the wave velocity and u Ni is the average 
velocity of the flat film. The interfacial shear stress is defined as 

T, = ^f s (U a+ V w )2, < 2 - 9 > 


where f s is interfacial fraction factor. 


Substituting equation (2.8) and (2.9) into equation (2.7) leads to 

)2 /-n ok 
(D - 2h) 3 


!*E* p 

dx p g‘s 


(Uan +u M) 2 ( D ~ 2h l^ 2 Si , 

t s — ( i + f ax >• 


( 2 . 10 ) 


The interfacial friction factor f s can be obtained from the pressure drop measurements and it is a 
function of the liquid flow rate. Equations (2.1)-(2.10) define the full equations of motion for the 
annular countercurrent gas-liquid two-dimensional falling film problem. We will start our 
analysis by considering the simplified equations for steady, laminar, one-dimensional smooth 
film motion with the countercurrent air flow , 
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u = 0 @ y = 0 , 


( 2 . 12 ) 
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9u 


^^ = T i @ y = h ’ 


(2.13) 


where tj represents the interfacial shear stress. 

The solution of equation (2.11) yields the following velocity distribution, 

4>h f f r, \ «2 ^ 


%(y) = 


\ + Jh-) y _yl 

P0hj y 


(2.14) 


where 


<t> = g- 


1 dp 
p dx 


(2.15) 


It follows from equation (2. 14) that the average velocity of the flat film is given by 


u Ni = 


^ , T i h 

3v 2fi 


(2.17) 


In order to simplify the two-dimensional model equations, we use a boundary layer type analysis. 
We scale the streamwise coordinate x by the unknown wavelength X, the normal coordinate y by 
h N ( Nusselt film thickness without any interfacial shear), the x-component of velocity by u N 
(average Nusselt velocity), the y-component of velocity by £un (where e = h N /X), pressure by 
pujM^ and time by Xu^. The scaled equations of motion and boundary conditions are 


5u 

ch 


da da dP 12 4 fd 2 u 2 d 2 u > ) 

+ u & +v ^=-^ + r + r\— +£1 —\’ 


(dy 2 


5x 2 


fdw dv dv\ 9P 4 

+ u “T~ + V -v . 1— - T= T — + FT 


dt^ u dK^ v dy -'^ + R 


'a * 2 ay 2 ' 


3u dv 


=o. 


(2.18) 

(2.19) 

( 2 . 20 ) 


y = 0 ; 


u = v = 0 , 


( 2 . 21 ) 
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y = h(x,t) ; 


v = h t +u h x . 


( 2 . 22 ) 


The tangential stress condition at y = h(x,t) becomes 

(ly + e2 ^)( 1 ' e2hx2 ) + 2 e2 (^' ^) hx = B 


(2.23) 


and the normal stress condition at y = h(x,t) is 


p p Be^ /3u jdv \ h x 8e 4 du h x 2 8e 2 dv 1 
• 8 + R ^ +E 3xJ 1+e 2 h 2- R 5T 1+E 2 h 2- R 5T 1+e 2 h 2 


+e 2 We 


A XX 


(l+e 2 h x 2 ) 3/2 


= 0, 


(2.24) 


where 


R = eRe, 


B = 


h N X i 

u nH 


(2.25) 

(2.26) 


We now consider simplification of the above model valid for large Reynolds numbers. 
For e « 1, Re » 1, eRe = R = 0(1) and We = 0(l/e) (valid for high surface tension fluids such 
as water) the model may be simplified by dropping the terms that are of order e 2 or higher. This 
leads to the extended form of the boundary layer model that is applicable in the presence of 
interfacial shear (in dimensional form): 


Du 1 dp 3 2 u 

DT = -pdx + g + V ^2 


(2.27) 


dp 

dy 


= 0, 


du dv 

5x + dy 


= 0 , 


(2.28) 

(2.29) 


u = v = 0 


@ y = 0, (2.30) 
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1! 

(2.31) 

d£_ 

dx 

- ah xxx +^ 

@ y = h, 

(2.32) 


v = ht + uh x 


@ y = h, (2.33) 


^P»_ nf 

dx P^s (D-2h) 3 


, , 2 5h . 

( 1 + fj 3x } ’ 


(2.34) 



(2.35) 


In this first approximation, the pressure is independent of y, the tangential stress condition is the 
same as that of the flat film and the normal stress condition includes surface tension and pressure 
gradient in the gas phase. In this paper we will only consider this boundary layer model. The 
second-order boundary layer model [8] will not be pursued because the pressure deviation across 
the film is negligible compared to the pressure applied on the interface because of the gas phase. 
The equations defining the boundary layer model are easier to deal with and expect this model to 
predict some observed wave characteristics as it includes all the leading order terms. 

3. Method of Solution 

We use the same formulation that can be found in Yu et al. [8]. to solve this problem, the 
only difference being in the tangential and normal stress boundary conditions. Because the flat 
film solution under the presence of tangential stress is different from that with zero shear stress, 
we need to find a different form of stream function that can deal with the flat film solution under 
the shear stress. We modify and write the stream function expansion as 



(3.1) 
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where b = . 

The first term of this stream function is the flat film solution under the presence of shear stress. 
The rest of the terms are used as correction terms to the base solution in the presence of waves. 
The no slip boundary condition at the wall can be automatically satisfied by this streamfunction. 
We use the same integral and boundary collocation method that was used in Yu et. al. On the free 
surface, there exists a general continuity condition for the tangential shear stress (written in the 
form of streamfunction), 

= b. (3.2) 

y=h 

Inserting the streamfunction expansion (1.3.1) into (1.3.2) results in the following equation, 

175h 2 a 5 + 120a 4 h + 63a 3 =0. 

In this equation, a 5 is linear with respect to the other variables so we can solve for a 5 as 

a 5 = -^2 (120a4h + 63a3) . (3.3) 

After substituting for a 5 into (1.3.1), we can rewrite the streamfunction in the following form, 

^(W) = a 2 (x,t)^l + ^)y 2 - ^ + a 3 (x,t)(y 3 - a 4 (x,t)^y 4 - ^ J 

- ^2 ( 120a 4 (x,t)h + 63a 3 (x,t)) (y 5 - ) . (3.4) 

The first condition used here is integral x-momentum equation to insure that the 
approximate velocity profile satisfies the global momentum balance for all times and locations. 
The x-momentum equation is integrated across the film, 




\0 


Y3u 9u 3u 1 dp 

S + u 5T + v 5T + 7di 



*-vH dy=0. 


(3.5) 


The result of integral x-momentum equation after simplification using Liebnitz' rule is 



Inserting the streamfunction expansion (1.3.4) into (1.3.6), the resulting equation may be 
normalized with respect to the flat film solution. This procedure gives the following equation, 
which was obtained from a symbolic manipulation program (MACSYMA), 

866250 ReWeH ^ + { 866250 ^ - [ 62300 A 4 2 H 6 + 221400 A 3 A 4 H 5 + 

( 148500 A 4 B + 297000A 2 A 4 + 188937 A 3 2 ) H 4 - ( 207900 A 3 B + 465300 A 2 A 3 ) H 3 
+ (144375 A 2 B+ 231000 A 2 2 )H 2 ] S - 

[ 288750 A 4 H 3 + 415800 A 3 H 2 + 288750 A 2 H]^- 

21 100 A 4 H 7 ^ - 46 140( A 3 ^ + A 4 ^) H 6 

3A a c 

- ( 41250 B + 82500 A 2 + 82500 A 4 + 102186 A 3 H 5 

3A-i ^Aj 3Ao 3Ao a 

- ( 86625 B + 115500 -3 f + 185625 A 2 + 185625 -^A 3 ) H 4 

3At 3At 3Ao - 3Ao „ 

- ( 144375 3X 2 B + 277200 + 346500 A 2 ) H 3 + 577500 H 2 } Re 

+10395000 H - 6930000 A 2 = 0. 


(3.7) 



In this equation, the normalized pressure in the gas phase and interfacial shear stress can be 
expressed as 




(U+ Ce) 2 ( 1- 2 a) 2 
(l-2aH) 3 


/i 2 5H , 

a(1+ r 5x 


(3.8) 


= RefV (U+Ce) 2 (1 -2a) 2 
8 p, s (1 - 2aH) 2 


where 


U_ 

U N 


a = 



,(a« 1) 



(3.10) 

(3.11) 

(3.12) 


For the problem considered here, D is 50.8 mm and h N is less than 0.5 mm when Re < 1500, so a 
is very small. For the range of values of H considered here(0 < H < 3), aH is less than 0.03. 
Using the approximation, 

— — — ~ 1+ na , lal « 1 , 
d-a) n 


we can simplify equations (3.8) and (3.9) to 



- m a (1 + 6aH)( 1 +j- ^ ), 

Dp 

B = ^m(l + 4aH), 


(3.13) 

(3.14) 


where 


m = ^L(U+Ce) 2 (l-2a) 2 . 
Pi 


(3.. 15) 



Taking the cross derivatives of the x and y momentum equations, equating like terms and 
evaluating near the wall leads to 


fa 3 ^ f a 4 v a 4v F ^ 

k ay 2 at v [ ay 2 ax 2+ 3y 4 JJ y=0 


(3.16) 


Inserting (3.4) into (3.8) results in the following equation (in dimensionless form), 


5Re 


3A 2 

3T 


3 2 A ? At 

' 40 ax 2 +144 H - 240A 4 = °- 


(3.17) 


The other collocation point we chose is at the interface, y = h. The x-momentum equation 
is evaluated at the interface (or infinitesimally close). At this position, the pressure gradient is 
specified by the capillary forces. The resulting partial differential equation, in term of the 
streamfunction, is 


a 2v F a 2v F a^ a 2 y ay 

3y3t + 3y3x c)y~ * 9y2 dx + pdx’^’ V 


a 3v p\ 



= 0. 


(3.18) 


Inserting the streamfunction expansion (3.3) into (3.18), the resulting equation was normalized 
with respect to the Nusselt flat film solution; it produces the following equation, 

1875 ReWeH |^| + ( 1875 ^ - [ 225 A 4 2 H 6 + 900 A 3 A 4 H 5 - ( 500 A 4 B + 1500 A 2 A 4 + 
864A 3 2 )H 4 -900( A 3 B + 3A 2 A 3 )H 3 - (1250 A 2 B + 1875 A 2 2 )H 2 )^ 

-( 1125 A 4 H 3 + 1800 A 3 H 2 + 1875 A 2 H)^- -75 A 4 H 7 ^- 

^ 3A, 3At f- 3A 4 3A A 3Ao 

180( A 3 W ‘ A 40X ^ H " ^ 1258 W + 375 A 2~SX + 375 a^T A 4 

+ 432 A 3 lx > « 5 - 3«X ^ B + | ^ + A 2 ^i + 3A 3 ^2) H< 

dA-y 3At dAy 1875 HR •> 

-[ 625-^ B + 900 -jf + 1875 A 2 B + 1875 A 2 ) H 3 

- 1875 H 2 } Re + 30000 A 4 H 2 + (27000 A 3 + 22500) H - 15000 A 2 = 0. 


(3.19) 



The continuity relation, 


du dv 

5T + 3y= 0 ’ 


is integrated across the film thickness, then using the Liebnitz' rule, the equation may be 
simplified and becomes 

■jL J ^ 

v(h) - v(0) - u(h) J udy = 0. (3.20) 

x o 


At the free interface, the velocities are related through the kinematic condition, 

V< h ) = U(h)| + |. 


(3.21) 


Substituting (3.20) into (3.21) and V(0) = 0 leads to 

I ♦£)«*- 


0 


ah , dT(h) n 
at + dx 


(3.22) 


When the streamfunction expansion (3.3) is inserted into (3.22), (3.22) becomes 

8 . 72 . ^ 4 . _ \dH 3H 


( o a m3 72 2 /t> 4 . . AdH dH 

^15 A 4 H -*t 5 A 3H +(B-t^ A 2)H J ax + ^ 




8 aA 


2dA' 


(3.23) 


The set of equations (3.7), (3.17), (3.19) and (3.23) can be solved for the unknown 
dimensionless variables A 2 , A 3 , A 4 and H as functions of position and time. 

The steady state solution of this set of four equations can be obtained by letting all the 
derivative be zero. This gives the algebraic equation, 

H 3 - — H 2 . i _ o 
n ss 2 n ss 1 ~ u ’ 


(3.24) 



where B is a constant and H ss is the dimensionless film thickness. Differentiating this equation 
with respect to H produces 
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H(3H - B) = 0, 

H = 0, orH = |. 

If B is negative, there is only one positive solution, if B is positive, there are two positive 
solutions. The solution of the cubic equation (3.24) for B < 0 is shown in Figure 2. From this 
figure, we see that, when the shear stress increases, the film thickness will also increase . 

4. Traveling Wave Model Formulation 

As discussed earlier, the wave reaches an asymptotic steady state after traveling several 
hundred film thicknesses. We are more concerned the fully developed waves than the linear 
growth wave region. So the traveling wave formulation will be used. 

Defining 


Z = X - Ce T, 


(4.1) 


where X, Ce, and T are as defined x/h N , c/u N and t/u N , motion with respect to the Z coordinate is 
steady, or 


_d d /0T\ a 
dz~[3z J T ax + (azJ T ar 


=o. 


(4.2) 


which allows the spatial and time derivatives to be written as 


_a_ __d_ 
ax ~dz 


(4.3) 


and 



(4.4) 





I £> 


or 

^O(H)-CeH)=0. (3.3) 

For the flat film, there is a relation, 

Y(H) - Ce H = 1 - Ce. (4.5) 

Substituting the streamfunction expansion (3.4) into above equation, the integral mass balance is 
used to solve directly for A 4 , 

A 4 = - (48A 3 H 3 +(75B+100A 2 )H 2 +150(Ce-CeH-l)). (3.4) 

Upon substitution of (7.3.4) into (3.7), (3.17) and (3.19) and transforming them to the 
traveling wave coordinate system produces the following equations: 

/dB dA 9 \ A f d 2 B d 2 A 9 ^ , , 3 9 

P + 2 l#J CeH4Re + \^l + 2 ^J H4 ‘ 288A 3 h3 -360(B+|A 2 )H 2 

+ 720 (Ce H - Ce + 1) = 0, (3.5) 


16 ReWeH 3 ^ + [16 H 3 ^ + (- B 2 H 4 - 2 B Ce H 3 + 12 (1 - Ce ) CeH + 36Ce 2 
72 Ce + 36) ^ - B ^|h 5 - 2 ^Ce H 4 + 6 ^|( Ce - 1) H ] Re+ 192(1 - 2A 3 )H 3 - 


( 960 B +1408 A 2 ) H 2 +1920 ( Ce H - Ce + 1) = 0, 


27720 H 3 Re We^- + {27720 H 3 - [180 A 3 2 H 6 + ( 576 A 3 B+ 1296 A 2 A 3 ) H 5 
- (3240 CeA 3 + 1 125 B + 2340 BA 2 + 2352 A 2 ) H 4 - ((4320 B + 2160 A 3 + 1 1040 A 2 ) 



Ce + 2160 A 3 ) H 3 + (- 12060 Ce 2 + (- 2160 B - 5520 A 2 ) Ce + 2160 B + 5520 A 2 ) H 2 

+ 39780 Ce 2 - 79560 Ce + 39780] ^ - 72 A 3 H 7 ^ - [ 144(A 3 ^| + B^ 3 )+ 

dAo dAo * dA 9 dB dA 0 

324( A 2 -^f + A 3"dz ^ H6 + n080Ce - ( 750 B + 780 A 2 ) - 780B 

dA ■) c dB dA -a dA a dA o * 

1568A 2 -^]H 5 +[(2160 gg - 1080 -^ + 5520 -^) Ce + 1080 -^] H 4 

dB dA 9 dB dA 9 ^ 

-[(2160^ + 5520 -g£) Ce + 2160 Jjg+5520-^ ] H 3 } Re 

+ 332640 H 3 - 221760A 2 H 2 = 0. (3.7) 


These three equations may be written as a first order system by using MACSYMA to simplify 
and rearrange the expressions. The resulting system is written in the general form: 

dH dOj 
dZ ~ dZ 2 ’ 

(3.8) 

d 2 H d0 2 
dZ 2 ~ dZ " °3 ’ 

(3.9) 

d 3 H d9o 

^ = - a f = f 1 (fi;Re,We,Ce), 

(3.10) 

CD 

A 

A 

(3.11) 

d 2 A 2 d0c 

dZ 2 = dZ =f 2 ( fi ;Re > We * Ce X 
dA. d0 fi 

= dZ = ^3^ ’ Re ’ ^ e ’ ^ ’ 

(3.12) 

(3.13) 


where the dynamic variable vector is defined by 



The functions, which are autonomous (no Z dependence) and highly nonlinear are given below 


fj(fi ; Re, We, Ce) = - ^ 3r — ((160 1 3 P z +0 2 (36Ce 2 +12Ce0 1 ( 1 - Ce) - 

20 1 3 BCe -72Ce-0 1 4 B 2 +36)+60 1 3 B z (Ce - l)-20 1 4 B z Ce-0 1 5 BB z )Re+192O GjCe- 
1 92OCe+0 1 2 (-960B- 14O80 4 )+0 1 3 (192-3840 6 )+192O), (3.15) 


tjffl ; Re, We, Ce) = - -j— 

loUj 


(0! 4 (B z +20 5 )CeRe+72O0 1 Ce-720Ce +80 / B^ 


+0 1 2 (-36OB-48O0 4 ) -2880 1 3 0 6 +72O), 


(3.16) 


f 3 (Q ; Re, We, Ce) = (2772O0 1 3 ReWef 1 +(2772O0 1 3 P z +0 2 (3978OCe 2 +0 1 2 (-12O6OCe 2 
+(-2 1 60B-5520 0 4 )Ce+2 1 60B +5520 0 4 >i- 0 1 3 ((432OB-2 1 60 0 6 + 1 1 0400^6 
+2 16006)+ 0^(3240 0 6 Ce- 1 125B 2 -2340 0 4 B-2352 0 4 2 )-7956OCe+0 1 s (-5760 6 B- 1296 
-18O0 1 6 0 6 2 +3978O)+0 1 3 ((-216O B z -552O0 5 )Ce+216O B Z +552O0 5 )+0 1 4 (216O B z 
+552O0(5))Ce+0 1 5 ((-75OB-78O 04)B Z -78O0 5 B-1568 0 4 0 5 ) 
+0 1 6 (-1440 6 B z -3240 5 0 6 )Re-22176O0 1 2 0 4 +33264O0 1 3 )/ 

((0 1 4 (lO8OCe-lO8O)lO8O0 1 5 Ce+0 1 6 (144B+324 0 4 )+72 0 1 7 0 6 )Re), (3.17) 

P z (Z)= - ma(l+6o(0 1 )( l+|-0 2 ), 

1 o 


(3.18) 



B(Z) = M (1 + 4a0 1 ) . 


(3.19) 
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Steady state solutions are determined by the setting all derivatives with respect to Z to 
zero, and solving the resulting equations. This produces the steady state vector, 

Q.ss = (eiss,o, o,e 4ss ,o,o) T , (3.20) 

where 

®lss = Hss, 

_ 3 

®4ss = ^2ss = 2^ss • 

A cubic polynomial for the steady state film thickness, H^, is obtained, 

h ss 3 + |h ss 2 - CeHss + Ce- !=°, (3.21) 


where 


Re 


B = m(l +4aH ss ) 


and 


M(1 +4 ocH^) 


= M (constant), 

\x Re 

M = -g-m. 


(3.22) 

(3.23) 


If B = 0 , then equation (3.21) is same as that studied by Yu et al. [8] and there are two solutions. 


and 


H (!) = 1 
n ss 1 



(3.24) 


(3.25) 



the first of which corresponds to flat film, and the second to another physical root existing for Ce 

3 g 

> ^ (which is from the transformation of the coordinate system). If B* 0 , we let D = y , and 

rewrite the cubic equation as follows, 

Hjg 3 + DH SS 2 - Ce Hjg + Ce - 1 = 0. (3.26) 

Depending on the values of Ce and D, this equation may have either zero, one, two or three 
feasible solutions for H $s . First we note that when Ce = 1, the roots are given by 

H ss = 0 ’ 

=^(-D + VD 2 + 4). (3.27) 


To determine the cusp point of equation (3.26) at which the three roots come together, we 
differentiate it with respect to H, once and twice and set the derivatives to zero. This gives 

3H ffi 2 + 20^ - Ce = 0 (3.28) 


and 


+ D = 0. 


(3.29) 


Solving the three simultaneous equations (3.26), (3.28) and (3.29) gives the cusp point 
coordinates. 


H SS = 






+ 1 = 3.1, 


Ce = 


Hss 3 + 2 
2 - 


-28.9 , 


D = 


(Hss~1) 2 (2H ss + 1) _ 
H ss (2-H ffi ) " 


The bifurcation set of equation (3.26) crossing which the number of solutions changes by two is 
obtained by solving equations (3.28) and (3.29) in a parametric form. 
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(3.30) 

(3.31) 

From these two expressions, we observe that 
D — > °° , Ce — ) 1 as H — > 0, 

D=0, Ce = 3 as H-» 1, 

D — ^ 00 , Ce — ) °° as H — > 2. 


Ce = ^ + 2 


2 - H 


ss 


D = 


(H,ss-1)^(2H SS + 1) 

H ss (2-H ss ) 


We show in Figure 3 the different regions in the (Ce, D) plane in which equation (3.26) has zero, 
one, two and three solutions. We now consider each of these regions: 

Region I: there are two film thicknesses for countercurrent or cocurrent gas flow (D positive or 
negative) 

Region II: In this region, there is no positive solution to equation (3.26). The reason may be 
when the gas flow is countercurrent, the film may be pushed away and eventually no film exists. 
Region III: In this region, Ce < 1 and there is only one positive film thickness for cocurrent or 
countercurrent gas flow. 

Region IV: In this region, there are three positive solutions if gas flow is the countercurrent and 
very high, possibly beyond the flooding transition where there is liquid flow in both directions. 



Ce = 1.0 


D = B/2 


Figure 3 The Bifurcation diagram of film thickness solution in the Ce and D plane. 


In this work, we consider countercurrent gas flow in region I, that is, B less than zero and 
Ce positive (before flooding), so there are two positive film thickness solutions. We plot the two 





solutions of the cubic equation (3.24) in Figure 4 for different values of B. When B = 0, the two 
solution branches intersect at Ce = 3.0, but if B < 0, the two branch will separate and form 
isolated branches shown in Figure 7.4. When the value of B is increased, the top branch will go 
up and the lower branch will go down. We now investigate the linear stability of these two 
separated solution branches. We represent the two steady states for this system as 


S-ss (1) =(Hss (1) .0, 0,§H a < 1 >,0,0), 

(3.32) 

e.ss (2) =(Hss (2) , 0,0, |H a ©, 0,0). 

(3.33) 


The characteristic equations that define the eigenvalues of these states are very lengthy 
and are not listed here. The polynomial equation was solved numerically for Re between 1 and 
1200, Ce between 1.0 to 3.2 for a fixed value of the Kapitza number, Ka (Ka = 3371, water). For 
each value of Ka, the unstable manifold dimensions for both steady states as functions of Re and 
Ce was generated. Figure 5 shows the unstable manifold dimension for the steady state at 
Ka = 3371. There is one Hopf line for the first steady state. The lower steady state (H s ®) is 
always unstable, that is, at least one real eigenvalue is positive. When Re is small, the Hopf line 
is very sensitive to the interfacial shear stress. When Re is increased, the celerity at the Hopf 
point is less sensitive. For the upper steady state film thickness H^ 1 ), the Hopf line drops 

quickly and reaches an asymptotic value that depends on the shear stress value for Re values of 
1200 or higher. For Ce greater than the Hopf celerity, the unstable manifold dimension is zero, 
meaning all small perturbations of the steady film are damped, and the film is linearly stable. 
Below the Hopf line, the unstable manifold dimension is 2, corresponding to a pair of 
eigenvalues crossing the imaginary axis. For the lower film thickness, the unstable manifold 
dimension is always great than one. 






The bifurcation diagram of the solutions shown in Figure 5 suggests all the interesting 
dynamic behavior occurs around the first flat film thickness. Additionally, the location of the 
wave inception line varies significantly with flowrate, agreeing with the experimental 
observations that fully developed films have different average wave velocities at different 
flowrates. The figure also shows that wave celerity is more sensitive to countercurrent gas flow 
for low liquid flowrates. The experimental data show that for liquid Reynolds number less than 
400, the wave celerity increases with increasing gas flow, but for Reynolds number larger than 
that the wave celerity changes little until the flooding point is approached [6]. 

5. Numerical Study of the Traveling Wave Model Equations 

The numerical method and analysis procedure used in this section are the same as those 
used in the last chapter. The aim of this study is to examine how the film profile character 
changes with varying wave celerities, Ce and shear stress B, for a fixed value of the Reynolds 
number and physical property group, Ka. 

The numerical study includes three parts: 

(1) determination of the capability of the model to produce statistically stationary wave 
forms and sensitivity to initial conditions, 

(2) a study of the character of the stationary wave forms, 

(3) a search for the values of Ce below which no bounded solutions exist. Integration of 
the system of ODEs required specification of initial conditions as well as a step size in the Z 
direction. Initial conditions were based on a perturbed flat film at the lower steady state 
thickness. In this study the film thickness was perturbed by 0.1% while all others retained their 
steady state values. 
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6. Sample Integration Results for Re = 600 and different values of interfacial shear stress 

In this study we are more concerned with the shear stress effect on the liquid film. 
Therefore, we use one typical flowrate corresponding Re = 600 and different values of shear 
stress. The dimensionless shear stress B is defined as 

B = M (1 + 4a H), (5.1) 


where 


M = 


mRe 


m = ^f.(U+Ce) 2 (1 -2 a) 2 . 

Pi 

The expression for the pressure deviation in the gas phase is 

dP g ,, . ov i 2 3H , 

rfy — ~ m cc (1 + 6aH)( 1 + ^ )• 


(5.2) 

(3.15) 


(3.13) 


Substituting (7.4.1) into (1.3.13) leads to 



8M , . 2 3H . 

' R^ a(1+6aH)(1+ f s dZ ) ‘ 


(5.3) 


In (7.4.3), the only unknown is interfacial shear stress f s which also appears in the expression for 
M. For the flowrate considered here, f s is about 0.02 [2], Because the value of B is not constant 
due to the wave present, a value of M is specified prior to each calculation. From expressions 
(4.1) and (4.3) the shear stress values in the normal and tangential direction are calculated and 
through expression (4.2) and (3.15) the dimensional gas velocity can be found to compare with 
the experimental data. For the flat film case, M is the same as B. 

A series of values of M are used to study the gas-liquid interaction at the interface. The 
values of M are in the range zero to -1.5 and the gas velocity associated with these values is from 
zero to 30 m/s. This range covers most laminar falling film flows. 

5.1 M = -0.2 



The celerity at the Hopf point for Re = 600 is 1.88. For Ce slightly below this value, the 
wave profile is sinusoidal with a single frequency. When the celerity is reduced, the wave 
amplitude become larger and the periodicity is lost. The transition from periodic to chaotic 
profile is similar to that with zero gas flow. We present here the wave profile for the smallest 
celerity value 1.785 for which a bounded solution exist. As shown in Figure 6, the wave 
thickness is larger than that for the free falling film due to the interfacial shear. The smallest 
value of wave celerity (1.785) for which a bounded solution exist has increased from that of the 
free falling film (1.73805). The film thickness and film slope phase plane shows the largest peak 
to substrate is 3.2 and the wave is more symmetric because the interfacial shear stress is in 
opposite direction to the gravity force(Figure 7). The power spectra and probability function 
shows no significant changes (Figure 8) other than the PDF is spread around the mean film 
thickness. This means the upward gas flow does not affect the nature of the falling film 
significantly at this value of M. 

.2 M = -1.0 

The celerity at the Hopf point for Re = 600 is 2.2. At near the Hopf celerity, the wave profile is 
sinusoidal with a single frequency. When reducing the celerity value the wave amplitude become 
larger and the also lose the periodicity. These transition is similar with the zero gas flow. The 
results presented here is for the smallest celerity value 1.91 that has bounded solutions. The wave 
thickness is larger than the free falling film due to the shear stress affect(Figure 9). The wave 
celerity 1.84 is also increased from 1.84 (M= -0.6) due to the increasing of the film thickness. 
The film thickness and film slope phase plane (Figure 10) shows the peak to substrate ratio is 
5.5 and the wave is more symmetry than before. The power spectra and probability function 
shows no significant changes (Figure 1 1). 

The celerity at the Hopf point for Re = 600 is 2.35. At near the Hopf celerity, the wave profile is 
sinusoidal with single frequency. When reducing the celerity value the waveamplitude become 
larger and the also lose the periodicity. These transition is similar with the zero gas flow. The 



results presented here is for the smallest celerity value 1.97 that has bounded solutions. The wave 
thickness is larger than the free falling film due to the shear stress affect (Figure 12). The wave 
celerity 1.97 is also increased from 1.91 (M= -1.0) due to the increasing of the film thickness. 
The film thickness and film slope phase plane (Figure 13) shows the peak to substrate ratio is 6 
and the wave is more symmetry than before. The power spectra shows no significant changes 
that means the gas flow does not affect the natural of the falling film flow and the center of 
probability function moves to the right direction (Figure 14). The PDF is also spreading along 
the film thickness and the gas flow increase the film randomness. 

Streamline maps are generated for different shear stress values to study the effect of 
countercurrent gas flow. For M = -0.6 there is no significant difference in the streamlines 
compared with that of the zero stress case [8]. There are two stagnation points on the interface 
streamline and there is a recirculation region under the wave peak (Figure 15). 

For M = -1.0, there is one recirculation region when the wave is not very large (Figure 
16(top)) . But a second recirculation region appears just under the interface for large waves 
(Figure 16(bottom)). This second recirculation is due to the interfacial shear stress and accounts 
for the increase in the interfacial mass transfer rate. 

For M = -1.3, there are two recirculation regions in the wave even with a peak to 
substrate ratio of 2.6 (top diagram of Figure 17). The size of the second recirculation region 
increases when the shear stress is increased while the size of the first recirculation region 
decreases (bottom diagram of Figure 17). 

7. Conclusions 

The models that developed for free falling film flow are extended to study the countercurrent 
gas-liquid falling film flow. Due to the complexity of the second order boundary model as well 
as the pressure variation across the film is small compared to the imposed gas phase pressure, the 
countercurrent gas flow affect was studied for the standard boundary layer model. A different 
stream fuction that can compencete the shear stress affect was developed and this stnean function 
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also can predict periodic solution. The descritized model equations were transformed to a 
traveling wave coordinate system. A stablity analysis of these sets of equations showed the 
presence of a Hopf bifurcation for certain values of the traveling wave velocity and the shear 
stress. The Hopf celerity was increased due to the countercurrent shear. For low flowrate the 
increases of celerity are more than for the high flow rate, that was also observed from the 
experiments. Numerical integration of a traveling wave simplification of the model also predicts 
the existence of chaotic large amplitude, non-periodic waves as observed in the experiments. The 
film thickness was increased by the shear. 
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Figure 7: Film slope vs. film thickness phase plane for Re = 600, M = -0.2 and Ce 
1.785. 
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Figure 9: Film thickness profile vs. distance for Re = 600, M = -1.0 and Ce = 1.91. 
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Figure 12: Film thickness profile for Re = 600, M = -1.3 and Ce = 1.97. 
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Figure 14: Normalized probability density function (top diagram) and power spectra 
(bottom diagram) for Re = 600, M = -1.3 and Ce = 1.97. 
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